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Finite Temperature Simulations from Quantum Field Dynamics? 

Mischa Salic, Jan Smit and Jeroen C. Vink* 

Institute for Theoretical Physics, University of Amsterdam 
Valckcnierstraat 65, 1018 XE Amsterdam, The Netherlands 

We describe a Hartree ensemble method to approximately solve the Heisenberg equations for the tp 4 model 
in 1 + 1 dimensions. We compute the energies and number densities of the quantum particles described by the 
tp field and find that the particles initially thermalize with a Bose-Einstein distribution for the particle density. 
Gradually, however, the distribution changes towards classical equipartition. Using suitable initial conditions 
quantum thermalization is achieved much faster than the onset of this undesirable equipartition. We also show 
how the numerical efficiency of our method can be significantly improved. 
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1. Hartree ensemble approximation 

Non-perturbativc equilibrium properties of 
quantum field theories are studied with great suc- 
cess using numerical lattice techniques. Similarly 
one would like to investigate quantum field dy- 
namics non-perturbatively. This has obvious ap- 
plications to e.g. heavy ion collisions and early 
universe physics. Furthermore, it could be used 
to study equilibrium properties in cases where 
Monte Carlo methods do not work, such as in 
QCD at finite density or chiral gauge theories. 
Using Monte Carlo techniques to compute the 
Minkowski time path integral is, of course, not 
possible. Neither can one directly solve the 
Heisenberg equations for the quantum field op- 
erators. Hence one has to resort to approxima- 
tions, such as classical dynamics || (for recent 
work see ||), large n or Hartree (see e.g. 
Here we discuss a different type of Hartree or 
gaussian approximation than used previously, in 
which we use an ensemble of gaussian wavefunc- 
tions to compute Green functions. Ref. jl) also 
contains a discussion of this method and a de- 
tailed presentation will appear elsewhere. 

We use the lattice tp 4 model in 1 + 1 dimensions 
as a test model, which has the following Heisen- 
berg operator equations, 



= 7T, 7T = (A — fl )0 — \<f> , 



(1) 



with A the lattice laplacian, \i the bare mass pa- 
rameter and A the coupling constant (we use lat- 
tice units a = 1). The Hartree approximation 
assumes that the wavefunction used to compute 
the Green functions is of gaussian form, such that 
the field operator can be written as 

(2) 



with a\ and a a time-independent creation and 
annihilation operators and time-dependent mean 
field tp and mode functions f a . All information 
is contained in the one- and two-point functions, 
which can be expressed (for gaussian pure states) 
as 



(<Px) = (fix, {<Px0y), 



rot rot- 
J x J y 



(3) 



The mode functions represent the width of the 
wave function, allowing for quantum fluctuations 
around the mean field. With Hartree dynamics 
we therefore expect to improve the classical dy- 
namics. Of course we cannot expect to capture 
all quantum effects, e.g. tunneling is beyond the 
scope of the gaussian approximation. 

The Heisenberg equations (|l|) provide self- 
consistent equations for the mean field tp and the 
mode functions f a , 
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In order to simulate a more general density ma- 
trix instead of the pure gaussian state assumed 
in eq. (|J), we average over a suitable ensemble of 
Hartrcc realizations by specifying different initial 
conditions and/or coarsening in time. In this way 
we compute Green functions with a non-gaussian 
density matrix p = J^iPiP?^ as 

i 

i 3 

It should be stressed that for typical Hartree 
realizations the mean field <p x is inhomogeneous 
in space. This is different from the ensemble av- 
erage of (p which may, of course, be homogeneous. 
We also note that the equations (Q) can in fact be 
derived from a hamiltonian. Since the equations 
are also strongly non-linear, this suggests that the 
system will evolve to an equilibrium distribution 
with equipartition of energy, as in classical statis- 
tical physics. 

2. Observables 

To assess the viability of our Hartree ensemble 
approximation, we solve the equations (^) start- 
ing from a number of initial conditions. We com- 
pute the particle number densities and energies 
from the connected two-point functions of ip and 
7r, after coarse-graining in space and time and av- 
eraging over initial conditions, 

^E^^Qconn = (6) 

-Y,e- ikz {Tr x+z 7r x ) coan = + (7) 

xz 

The over-bar indicates averaging over some time- 
interval and initial conditions as in eq. (||), tu^ is 
the energy, the (number) density of particles 
with momentum k and N the number of lattice 
sites. 

For weak couplings, such as we will use in our 
numerical work, the particle densities should have 
a Bose-Einstein (BE) distribution, 

n k = l/(e Wfc / T -l), u>l = ?7i 2 (T) + 2-2cos(fc),(8) 
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Figure 1 . Particle number densities as a function 
of u> at various times. The straight line gives a 
BE distribution with temperature T/m = 1.7. 
(The lattice volume Lm = 32, coupling X/m 2 = 
0.083 and inverse lattice spacing 1/am = 8; m = 
m(T = 0)). 

with T the temperature and m(T) an effective fi- 
nite temperature mass. We have verified that this 
is indeed the case using Monte Carlo simulations 
of the cuclidian time version of the model. 

3. Numerical results 

First we use initial conditions which correspond 
to fields far out of equilibrium: gaussians with 
mean fields that have just a few low momenta 
modes, 

jmax 

(p x = v, n x = A "E cos(kjX + ctj). (9) 
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The ctj are random phases and the A is a suitable 
amplitude. Initially the mode functions are plane 
waves, e lkx /v2ujL, such that all energy resides in 
the mean field. 

The results in Fig. |l| show that fast, well be- 
fore tm « 200, a BE distribution is established 
for particles with low momenta. Slowly this ther- 
malization progresses to particles with higher mo- 
menta, while the temperature T k, 1.7m remains 
roughly constant. Such a thcrmalization does not 
happen when using a single Hartrcc realization 
with a homogeneous mean field ip. The differ- 
ence may be understood, because particles in our 
method can scatter via the inhomogeneities in the 
mean field. 



3 



Log (1+1/n) 




Figure 2. Particle number densities as a function 
of ui, at early times (a and b) and large time (c). 
(Lm = 25.6, A/to 2 = 0.5 and l/am = 10). 

Next we want to speed-up the thermalization of 
the high momentum modes. Therefore we use dif- 
ferent initial conditions in which the energy is dis- 
tributed more realistically over the Fourier modes 
of the mean field. The modes are the same as be- 
fore but mean fields arc drawn from an ensemble 
with a BE-likc probability distribution, 

P (ip k ,n k )^ C M~(^ k/To ~l)(^l+tolipl)/2co k ].(10) 

We also want to probe large time scales, which 
were not yet reached in the weak coupling simu- 
lation of fig. [j]. 

Fig. ||a shows the particle density computed 
without the mean field part in eq. (Q). Already 
at tni = 20 a credible BE distribution has estab- 
lished, with a temperature that increases until it 
saturates at T « 1.0m at tm = 100 (T = at t = 
0). Fig. ^) shows the densities at tm = 100, now 
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Figure 3. Particle number densities at large 
times, tm = 25600 and high initial temperature, 
Tq = 5m. (Lm = 25.6, X/m 2 = 0.5 and l/am = 
10). 

with the mean field contribution included. The 
straight line is a BE distribution at T w 1.01m. 

At very large times, Fig. ||c, the BE distribu- 
tion persists in the large 10 region, albeit with a 
slightly higher temperature T w 1.14m, but for 
small u> one clearly sees deviations. The parti- 
cle densities here can be fitted rather well with 
the classical form n — T c i/lj (dashed curve). 
The corresponding classical temperature slowly 
decreases, as energy disperses more and more to 
the large momentum modes. It is remarkable that 
the "decay time" of this temperature is some two 
orders of magnitude larger than the warming-up 
time shown in Fig. ^|a. 

Finally, in Fig. we show data for a large time 
tm w 25000 and at a higher (initial) temperature, 
To = 5m. Here cquipartition has extended to 
all modes with w^12to, as is evident from the 
accuracy of the fitted curve n = 1.9m/w (dashed 
line). 

Even though the timescales of relaxation to 
a BE distribution and relaxation to classical 
cquipartition appear to be widely separated, the 
simulations of Fig. || could not easily be used 
to study finite temperature equilibrium physics: 
only around tm = 100 — 200 there may be a win- 
dow in which the fields have thermalizcd to a BE 
distribution, without significant deviations due to 
the up coming classical equipartition. To increase 
such a window we did a simulation at weaker cou- 
pling and higher temperature. Now we find den- 
sity distributions as shown in Fig. 0. Already af- 
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Figure 4. Particle number densities as a func- 
tion of a;, starting from BE type initial conditions. 
The straight line through the origin is a BE distri- 
bution with temperature T/m = 6.7. (Lm = 23, 
A/to 2 = 0.083 and 1/am = 22). 

ter a relatively short time, before tm « 100, the 
particles have acquired a BE distribution up to 
large momenta uj/m w 12. Note that even with 
the BE-type initial conditions, the fields initially 
are out of equilibrium: energy is initially carried 
by the mean field only but within a time span of 
tm « 200 it is redistributed over the modes. 

Finally we try to improve the efficiency of our 
method. Since we have to solve for N mode 
functions, on a lattice with TV sites, the CPU 
time for one time-step grows oc ./V 2 . However, 
mode functions corresponding to particles with 
momenta much larger than T should be irrele- 
vant, since these particle densities are exponen- 
tially suppressed. This suggests that we can dis- 
card such modes. This is tested in Fig. ||, where 
we compare simulation results using all mode 
functions with results obtained using only a quar- 
ter of the mode functions. This corresponds with 
mode functions with initial plane wave energy 
uj^Hm. Clearly the results for particles with sig- 
nificant densities n are indistinguishable. For par- 
ticles with energy larger than f=a 17m, there are no 
longer mode functions that can provide the vac- 
uum fluctuations and consequently the particle 
density defined by (0) drops to —1/2. 

4. Conclusion 




Figure 5. Hartree dynamics with all mode func- 
tions (drawn lines) and with a reduced number 
(32 out of 128) mode functions (dots). (Lm = 5.7, 
A/to 2 = 0.083 and 1/am = 22). 

quantum thermalization in a weakly coupled 
scalar field model in real time. Only after times 
much longer than typical equilibration times, the 
approximate nature of the dynamics shows up in 
deviations from the BE distribution towards clas- 
sical equipartition. We have furthermore found 
that these simulations can be done efficiently us- 
ing a limited number of mode functions. 
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Using a Hartree ensemble method, we have 
demonstrated that we can simulate approximate 



